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Abstract: Production of muons and neutrinos in cosmic ray interactions with the 
atmosphere has been investigated with Monte Carlo models for hadronic interactions. 
The resulting conventional muon and neutrino fluxes (from n and K decays) agree 
well with earlier calculations, whereas our prompt fluxes from charm decays are 
significantly lower than earlier estimates. Charm production is mainly considered as 
\ a well defined perturbative QCD process, but we also investigate a hypothetical non- 

perturbative intrinsic charm component in the proton. The lower charm rate implies 
better prospects for detecting very high energy neutrinos from cosmic sources. 



^ ■ 1 Introduction 

<D ' 

■ The flux of muons and neutrinos at the earth has an important contribution from decays of 

particles produced through the interaction of cosmic rays in the atmosphere (for a recent intro- 
duction see ]l]]). This has an interest in its own right, since it reflects primary interactions at 
energies that can by far exceed the highest available accelerator energies. It is also a background 



in studies of neutrinos from cosmic sources as attempted in large neutrino telescopes, such as 
Amanda Q, Baikal ||, Dumand j| and Nestor ||. 

Here we present a detailed study of muon and neutrino production in cosmic ray interactions 
with nuclei in the atmosphere using Monte Carlo simulations ||. 

At GeV energies the atmospheric muon and neutrino fluxes are dominated by 'conventional' 
sources, i.e. decays of relatively long-lived particles such as tt and K mesons. This is well 
understood from earlier studies ]7|, ||, [|, with which our investigations agree. With increasing 
energy, the probability increases that such particles interact in the atmosphere before decaying. 
This implies that even a small fraction of short-lived particles can give the dominant contribution 
to high energy muon and neutrino fluxes. These 'prompt' muons and neutrinos arise through 
semi-leptonic decays of hadrons containing heavy quarks, most notably charm. 

Available data in the multi-TeV energy range, obtained with surface and underground detec- 
tors (see e.g. refs. [10-13]), are still too discrepant to draw definitive conclusions on the flux of 
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prompt muons and neutrinos from charm. Furthermore, estimates of these prompt fluxes [7,14— 
21] vary by a few orders of magnitude due to the different models used to calculate the charm 
hadron cross section and energy spectra. This huge model dependence is due to the need of ex- 
trapolating charm production data obtained at accelerator energies to the orders-of-magnitude 
higher energies of the relevant cosmic ray collisions. Obviously, this extrapolation can only be 
trustworthy if starting from proper charm production data and using a sound physical model. 
The main contribution of our study is in this context. First, we use recent charm production 
cross section measurements that form a consistent set of data, but disagree with some of the 
early measurements that were substantially higher. Secondly, we apply state-of-the-art mod- 
els to simulate charm particle production through perturbative QCD processes in high energy 
hadron-hadron interactions. In addition, we investigate a possible non-perturbative mechanism 
using the hypothesis of an intrinsic charm quark component in the nucleon. 

In the following we first (section 2) discuss the generalities of cosmic ray interactions in the 
atmosphere resulting in a set of transport or cascade equations for particle propagation. These 
equations are then solved by two different methods: a direct Monte Carlo simulation of the 
cascade interactions (section 3) and a semi-analytic method (section 4) giving consistent results 
for the conventional and prompt muon and neutrino fluxes. Section 5 gives an account of the 
Monte Carlo model used to obtain the energy spectra of secondaries in the basic hadron-hadron 
interaction, in particular concerning charm production in perturbative QCD. In section 6 we 
investigate consequences of the non-standard hypothesis of a non-perturbative intrinsic charm 
quark component in the nucleon. We then (section 7) compare our results with previous model 
calculations and discuss differences in terms of the different charm production models. We 
conclude (section 8) by some remarks and by putting our results in a general context of various 
astrophysical sources of high energy neutrinos. 



2 Cosmic ray interactions in the atmosphere 



2.1 The spectrum of cosmic rays 

Fluxes of secondary particles (hadrons and leptons) originate from nucleon-nucleon encounters, 
even when the nucleons are bound in nuclei, because nuclear binding energies are much lower 
than the energies of interest in this study (100 GeV - 10 9 GeV). So the relevant quantity to 
consider is the flux of nucleons. Following Jl|, [| 17 1 we have assumed a power law primary 
nucleon flux 



->JV 



(E) 



nucleons 



cm ssr GeV/ A 



1.7 E~ 2 - 7 for E < 5 • 10 6 GeV 
174 E~ 3 for E > 5 • 10 6 GeV. 



(1) 



The normalisation constant 1.7 is derived [22] from the directly measured primary spectrum 
using balloon-borne emulsion chambers in JACEE p3[. To within some 10% this agrees with 



more indirectly derived spectra based on measured atmospheric muon fluxes [24|, and is also 
compatible with the data discussed in ref. |25| ]. The cosmic ray composition is dominated by 
protons with only a smaller component of neutrons in nuclei [^, . Only primary protons are 
considered here, since in this study we are interested in quantities that are essentially indepen- 
dent of the cosmic ray composition. At the energies of interest (-E^lOOGeV), the cosmic ray 
flux can be considered isotropic (the anisotropy being ^5% pG|]). 
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2.2 The model for the atmosphere 



In studying the propagation of particles through the atmosphere, an important quantity is the 
amount of atmosphere X, in g/cm 2 , traversed by the particle. This so-called slant depth is 
the integral of the atmospheric density from the top of the atmosphere downward along the 
trajectory of the incident particle. At distance i from the ground along a direction at an angle 
9 from the zenith, the slant depth is defined as 

roc 

X(£,9) = / p[h(l,e)]dl, (2) 
h 

where p[h(l, 9)] is the atmospheric density at the altitude h(l,9), 

i / 2 

h{l, 9) = JRl + 21 R® cos 9 + P - R® « I cos 9 + — — sin 2 9. (3) 

Here R® is the radius of the Earth and the approximate equality applies for zenith angles not 
too far from vertical, 0^60°. 

The atmospheric density is not a simple function of height, even neglecting local atmospheric 
turbulence. The temperature, which is related to the density through the equation of state, de- 
creases with increasing height until the tropopause (8-17 km), stays almost constant in the lower 
stratosphere (— 56.42°C up to 20-30 km), then increases until the stratopause (50km) before de- 
creasing again at the highest altitudes (> 50 km). However, since most particle interactions 
occur at heights between 10 and 40km (demonstrated in Fig. p] below), we need only a simple 
model for the density profile of the stratosphere. We therefore adopt an isothermal model, 

p(h) = P oe- h ' h \ (4) 

with scale height Hq = 6.4 km and Xq = poh$ = 1300 g/cm 2 , values which adequately describe 
the density of the stratosphere (< 2% error in the vertical depth between 10 and 30 km and 
< 16% between 30 and 40 km). 

Concerning the atmospheric composition, a good approximation, valid up to a height of 
100km, is 78.4% nitrogen, 21.1% oxygen and 0.5% argon (obtained from data in (27]]). This 
leads to an average atomic number of (A) = 14.5. 



2.3 Particle interactions with air nuclei 

To obtain the flux of atmospheric muons and neutrinos one needs to consider the particle pro- 
duction mechanisms in strong interaction dynamics. The cosmic ray particles, represented by 
protons (see section 2.1), interact with nuclei in the atmosphere to produce secondary parti- 
cles. These proton-nucleus collisions can, for our purposes, be well represented by the simpler 
proton-nucleon collisions and a rescaling of the cross section 

a{ P A) = A a a(pN) (5) 

using a power dependence on the number A of nucleons in the target nucleus. For inclusive 
cross sections, with a{pN) of order 10 mb, the interaction occurs with nucleons at the surface 
resulting in a ~ 2/3 as verified experimentally. 

The inelastic pN interaction produces secondary hadrons with a multiplicity increasing es- 
sentially logarithmically with the cms energy (s p n = 2ra^(?E v ). The formation time of a hadron 
is the normal strong interaction time scale. In the particles rest system this corresponds to a 
formation length of ~ 1 fm, which is Lorentz transformed with a 7-factor to the target nucleus 
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rest frame and thus becomes proportional to the energy of the particle |28|| . Therefore, fast 
particles have formation lengths that exceed the size of the nucleus whereas only slow particles 
are formed and can re-interact within the target nucleus. Therefore, intra-nuclear cascade effects 
are not important for the energetic particle production studied here. 

Of importance for our considerations are the energy distributions of secondary hadrons pro- 
duced in the collisions 

dn(kA^hY;E k ,E h ) _ 1 dajkA hY; E k , E h ) 

dE h a kA (E k ) dE h U 

where dn(kA — » hY; E k , E^) is the number of hadrons h with energies between E^ and E^ + dE^ 
produced in the collision of the incoming particle k with an air nucleus of atomic number A, and 
a k A is the total inelastic cross section for particle k - nucleus A collisions. Experiments studying 



proton-nucleus [29| and heavy ion |3(J collisions obtain energy spectra that are approximately 
the same as in proton-proton collisions, confirming that the interactions are essentially proton- 
nucleon. So we adopt an energy distribution dn k h/dEh independent of the atomic number of 
the target. 

The energy spectra in Eq. @ can also be expressed as 

dn kh _ 1 da(kA -> hY) 
dxF a k A dxF 

in terms of the scaled longitudinal momentum or Feynman-x variable xp = p z /Pz,max (~ Eh/E k 
at large energies) where the z-axis is along the incoming particle momentum. If these distri- 
butions are independent of the cms energy (i.e. incoming particle energy Ek), then 'Feynman 
scaling' holds. The validity or breaking of this scaling in different models for particle production 
is an important issue as will be demonstrated later. 

To obtain the energy spectra of the particles produced in proton-nucleon collisions we use 
the Lund Monte Carlo simulation programs Pythia and Jetset [pi]] . These have proven very 
successful in describing the multi-particle final state in various kinds of interactions, including 
hadron-hadron collisions. An advantage with this Monte Carlo approach is the access to the 
complete final state as well as a proper account of the decay of unstable particles. Conventional 
muons and neutrinos are obtained from an inclusive event sample generated with Pythia in 
a mode simulating minimum bias proton-proton interactions (including diffractive scattering). 
The particle production results from Lund model [32] hadronization of colour string fields be- 
tween partons scattered in semi-soft QCD interactions. The prompt muons and neutrinos, on 
the other hand, are obtained from a dedicated charm production simulation using Pythia. 
Here, charm particles arise from the hadronization of charm quarks produced in the processes 
gg — > cc and qq — > cc as calculated with leading order perturbative QCD matrix elements. A 
more detailed account of the Monte Carlo model is given in section 5. Since non-perturbative 
charm production is neither well established nor well defined, it is not part of our main Monte 
Carlo study but investigated separately based on the intrinsic charm hypothesis in section y. 



2.4 Particle propagation in the atmosphere 

Propagation of high energy particles through the atmosphere may be described by a set of 
transport or cascade equations. In principle, the transport equations for nucleons, mesons, 
unstable baryons and leptons are coupled, but under the reasonable assumptions made below 
they can be greatly simplified. 

Nucleons constitute the initial primary flux. We consider nucleon absorption and regenera- 
tion in nucleon-air inelastic collisions, but neglect the certainly small contribution to the nucleon 
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flux from the interaction of unstable hadrons with air nuclei. Absorption is described by the 
interaction thickness A at of nucleons ./V in air, i.e. the average amount of atmosphere (in g/cm 2 ) 
traversed between successive collisions with air nuclei. It is given by 

X N (E) = — ^y-Jr — , (8) 

T,a VNAiE) n A (h) 

where n^(/i) is the number density of air nuclei of atomic number A at height h and ctna{E) is 
the inclusive inelastic cross section for collisions of nucleons with nuclei A. Note that to a good 
approximation \^{E) does not depend on the height h because the atmospheric composition is 
approximately independent of the height up to 100 km. 

Nucleon fluxes develop according to the cascade equation 

^» = -pL + S (NA^NY). (9) 
dX Ajv 

Here (j)]\r(E, X, 0) is the nucleon flux at slant depth X in the atmosphere at zenith angle 0, Xj^(E) 
has been defined in Eq. (||), and S(NA — > NY) is the nucleon-nucleon regeneration function in 
air 

S(NA NY) = r dtf dn(NA —> NY:E',E) 

Je an{E') dE 

Mesons and unstable baryons, in addition to interact with the atmosphere, can also decay. 
The decay length d(E), i.e. the distance traveled in a mean decay time, is simply 

d(E) = c/3jT (11) 

where r is the particle proper lifetime, 7 = E/mc 2 is its Lorentz factor, m its mass, and (3 its 
speed in units of the speed of light c. The decay length increases with particle energy because 
of relativistic time dilation; faster particles can travel longer before decaying. This implies an 
increased probability to interact before decaying. It is exactly because of this energy-dependent 
competition between decay and interaction that the muon and neutrino fluxes from charm 
mesons overcome those from pions and kaons at high enough energy. 

We assume that mesons and unstable baryons (collectively unstable hadrons) are generated 
in nucleon-air collisions and regenerated in hadron-air collisions, but neglect generation of 
unstable hadrons of other types in collisions of hadrons against air nuclei. This approximation 
is reasonable since the fluxes of unstable hadrons are at least a factor of ~ 10 smaller than the 
fluxes of nucleons. Thus for mesons and unstable baryons we have 

% = S(NA ^MY)-^-^ + S{MA - MY), (12) 
dX pdM am 

where Xm(E) is the hadron interaction thickness in air (analogous to Eq. (8)), and S^iV^l — > 
MY) and S(MA -> MY) are defined analogously to Eq. @. 

Finally we consider muons, muon-neutrinos and electron-neutrinos. At the energies we are 
interested in, energy loss, absorption and muon decay can be neglected and the transport equa- 
tion for lepton £ = /i^, u^, v^, u e , v e contains only source terms 

jH=5XM-*n (13) 

M 
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where the sum runs over all mesons and unstable baryons decaying into muons and neutrinos 
M = tt ± , K±, K°, D ± , D°, D°, Df, Af. S(M -» IY) describes lepton production in hadron 
decays, 



S(M -+ tY) = r dE M MM ~* ^ E »> ^ 

JE pd M {E M ) 



dE 



(14) 



where dn(M — > £Y; Em , E) is the number of leptons with energy between E and E + dE 
produced in the decay of hadron M. 

Muons and neutrinos born out of pions and kaons are traditionally called 'conventional,' 
while those born out of charmed hadrons are called 'prompt.' This originates from the fact that 
up to ~ 10 7 GeV very short-lived charmed particles have negligible probability of being absorbed 
in the atmosphere before decaying. Up to ~ 10 7 GeV the prompt flux is therefore essentially 
independent of the zenith angle. The conventional muon and neutrino fluxes are instead lower in 
the vertical direction, where the amount of atmosphere traversed in a given meson decay length 
is larger. The prompt flux is therefore relatively more important in the vertical direction, and 
we will predominantly consider this direction. 

3 Simulation of cascade interactions 

One way of solving the transport equations described in the previous section is to simulate the 
particle cascade with a Monte Carlo program. Here we describe our simulation algorithm. 

A cosmic ray proton is generated. Its energy is drawn from a flat distribution in \ogE, and 
a weight is assigned to it in order to reproduce the shape of the primary spectrum. 

An interaction height h for the cosmic ray proton is then chosen in the following way. Primary 
nucleons propagate down through the atmosphere according to Eq. @ without the regeneration 
term S(NA — > NY). From the solution to this equation, (f)(h) = <f>ao exp (— X(K)/\n), the 
probability distribution for the primary interaction height can be obtained. Using a standard 
Monte Carlo technique, we generate this distribution by replacing </>(/i)/^oo with a uniform 
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Figure 1: Distribution of the altitude for the primary interactions as obtained in the cascac 
simulations. 
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random number R G]0, 1[ and then solving for the interaction height h. This can be done 
analytically for the isothermal atmospheric model in Eq. (4) in the vertical direction, for which 
we obtain 

h = -ho In (^^) • (15) 

Fig. [l] shows the height distribution so obtained (neglecting the logarithmic energy dependence 
of A at), which confirms that, under the assumptions made, most particle interactions occur at 
heights between 10 and 40 km. 

A proton-nucleon interaction is then generated in full detail with Pythia |3lJ] resulting in a 
complete final state of particles. Secondary particles are followed through the atmosphere where 
they decay or interact producing cascades. Secondary nucleons give a flux that is rather small 
compared to the primary flux and could therefore be neglected as a first approximation. To be 
more precise, we do include the main effect of this correction by taking into account secondary 
nucleons that have an energy of at least 30% of the primary one. Nucleons with a lower energy 
give a negligible contribution compared to the primary flux due to its steep energy spectrum 
Eq. (|H). These leading nucleons emerging in the interactions are therefore allowed to generate a 
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Figure 2: The E 3 -weighted flux of muons + pT ), muon-neutrinos (v„ + v^) and electron- 
neutrinos (u e + v e ) from decays of the specified particles. The error bars indicate the statistical 
precision of the Monte Carlo simulation. 
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secondary interaction at a height 



-h ln(e- H ^_h^L) } (16) 



X 

obtained analogously to Eq. (|l5|) but taking into account the finite height H of the primary 
interaction. The procedure is iterated until the energy of the leading nucleon from an interaction 
falls below 30% of the primary cosmic ray proton energy. 

Secondary mesons and unstable baryons are traced through the atmosphere until they either 
decay or interact. Which of these occurs is decided by comparing simulated decay and interaction 
lengths 

L dec = -d M {E) In R x (17) 

and 

L tnt = H + h ln{e- H ^-^-^y (18) 

where R\ and i?2 denote uniform random numbers G]0, 1[ and H is the height at which the 
traced particle has been produced. The decay length cIm(E) and the interaction thickness Am 



are given in Eqs. (11) and (|§|) respectively, the atmospheric scale height ho and depth Xq are 



defined in sect. 2.2. Eq. fliq ) is obtained in a way analogous to Eq. (16). 

Particle decays are fully simulated with daughter particle momenta. In case of interactions, 
the interacting particle is regenerated in the same direction but with degraded energy, chosen 
according to the appropriate leading particle spectrum. Considering only the most energetic 
'leading' particles in secondary interactions is justified because they give the dominant contri- 
bution to the lepton fluxes. Moreover, other particles with lower energy are much fewer than 
the particles of the same type and energy produced in primary interactions. 

The particle decay-interaction chain is then repeated until all particles have decayed, have 
hit the ground or their energy has fallen below the minimum energy of interest, 100 GeV. Energy 
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Figure 3: The E 3 -weighted vertical flux of muons, muon-neutrinos and electron-neutrinos from 
conventional (ir,K decays) and prompt (charm decays) sources and their sum ('total'). The 
solid lines are from the cascade simulation (section 3) and the dashed lines are from the analytic 
Z -moment method (section 4)- 
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No 7 A E i % 


Conventional fj,~ + /i + 
Prompt [i~ + /i + 


2.0 • 1CT 1 1.74 7.0 • 10" 3 5.3 • 10 5 2.10 2.2 • 10 1 
1.4 • 10~ 5 1.77 2.8 • 10~ 8 9.2 • 10 5 2.01 4.3 • 10~ 4 


Conventional v„ + Vn 
Prompt z/ M + Ufj, 


1.2 -10 -2 1.74 2.0 -lO -3 6.3 • 10 5 2.17 3.8 • 10 u 
1.5 • 10" 5 1.77 3.1 -10" 8 1.2 -10 6 1.99 3.1 • 10~ 4 


Conventional v e + v e 
Prompt v e + u e 


4.2 ■ 10" 4 1.63 7.0 • 10~ 3 5.0 • 10 4 2.12 8.4 • 10~ 2 
1.5 • 10" 5 1.77 3.0 -10" 8 1.2 -10 6 2.02 4.9 • 10" 4 



Table 1: Values of parameters in Eq. 
cascade simulations in Fig.^. 



lW obtained from fits to the Monte Carlo results of the 



spectra for muons and neutrinos are finally obtained by counting muons and neutrinos with the 
initially- assigned primary proton weight. 

The resulting fluxes of muons and neutrinos from different parent particles are shown in 
Fig. [2]. For charmed particles the figure clearly demonstrates the dominance of the D ±,a mesons, 
while for conventional fluxes the dominant source varies with the type of lepton considered. 
Summing the various contributions gives the inclusive fluxes in Fig.|3|. As can be seen, the 
prompt contribution from charmed particles dominates at high energies. 

The results for the inclusive prompt and conventional fluxes can be parametrised as (similarly 
to @ , cf. with the primary flux and Eq. (j3g|) below) 

N E-^- 1 /(l + AE), E<E , 

N' E-^- x /{l + AE), E>E , 
with an accuracy of typically better than 10% using the fitted parameter values in Table |l[ 



HE) 



(19) 



4 Approximate analytic solutions 

Approximate analytic expressions for the muon and neutrino fluxes can be found from the cas- 



cade equations in sect. |2.4| by interpolation of high-energy and low-energy asymptotic solutions. 
This is done in the standard treatment for power law primary spectra and scale-invariant interac- 
tion cross sections jl], 0, ||] . We wish to generalize the standard treatment to include non-scaling 
effects. 

The cascade equations for nucleons and mesons (and unstable baryons) can, using Eqs. (9,10, 
12), be written 



3JV 



dX 

d^M 
dX 



+ Zmn\ - 
An an 



pd 



M 



<t>M . ry <t>M . 

t r 6MM\ V ^NM\ — , 

Am Am Ajv 



(20) 



(21) 



where the spectrum- weighted moments for generation Znm and regeneration Znn, %MM in 
hadronic collisions are generally defined as 



■>kh 



dE' 



: (E',X,6) X k (E) dn(kA — > hY; E' , E) 



h(E,X,6) X k (E>) 



dE 



(22) 



It is assumed that the fluxes of nucleons, mesons, unstable baryons, muons and neutrinos can 
be approximated in the factorized form <fii(E, X,8) = E~@' <f>i(X,9), with appropriate values of 
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the exponents 0i in the low- and high-energy asymptotic regimes. We consider both a primary 
spectrum (j)]y(E) oc E~^~ l with a constant spectral index 7, and the primary spectrum with 
a knee as given in Eq. (|l|). Since nucleon and meson fluxes develop rapidly in the atmosphere, 
their ratios are essentially independent on the depth X. To a good approximation one therefore 
obtains the energy- dependent spectrum-weighted moments 

Z kh (E) - J e dE {-) — , (23) 

which we will estimate numerically. In previous studies it has been assumed that Feynman 
scaling holds, such that the distributions in xp = E^jE^ of produced particles are energy 
independent (cf. Eq. (|7|)). With this additional assumption the Z-moments become simply 

z sealin 9 = j\ F ^dx F . (24) 

We will not make this scaling assumption, but investigate its validity by our calculation of the 
energy-dependent Z-moments. 

An approximate solution for the nucleon fluxes is then 

<t>N = e~ X/AN <j> N (E), (25) 
where the nucleon attenuation length A^r is defined as 

An(E) = X f E \ (26) 

and <Pn{E) is the primary nucleon spectrum. 

Concerning mesons and unstable baryons, at sufficiently low energies the interaction term 
can be neglected and so also the regeneration term. One then finds 

^ = E^. e -X/As ME) . (27) 

1 — Ann Ajv 

In the high energy regime it is the decay term that can be neglected, and one finds in a similar 

way 

u- v, Zmm e ~ x / A M _ e -x/A N 

^ i-zV 1 - Aat/Aju ^ < 28 > 

with the meson attenuation length Am defined analogously to Eq. (p6[). Notice that at high 
energies the spectral index of the meson flux is the same as that of the primary flux, <f>^? oc 
while at low energies the meson spectrum is flatter by one power of energy, (jyffi oc E^ 1 , 
because of the implicit proportionality of the decay length dM to the energy E. 

The spectrum-weighted moments for hadron generation in hadron-air collisions can be 
rewritten as 

*>m = r**^ (W 1 in{kA ^ E «- E *\ (29, 

JE h cr kA (E h ) \E h J dE h 

The numerical evaluation of these Z-moments was made by applying the prefactor in the inte- 
grand to the hadron spectra dn^/dE^ generated at different incoming energies E k between 10 2 
and 10 9 GeV with the Pythia Monte Carlo (see section 5 for details on the generation mecha- 



nism). Total inelastic cross sections a k A{E) in Eq. (|29| ) were taken from ref. [33] when available. 



The spectra of regenerated kaons and D mesons, which cannot be used as beam particles in 
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Figure 4: Energy- dependence of production Z '^-moments, Eq. for incoming particle k pro- 
ducing hadron h. Solid lines are the results of our model using the initial spectrum with a 'knee', 
Eq. whereas the dotted lines are obtained with a constant spectral index 7 = 1.7. Dashed 
lines show the values of j^J based on Feynman scaling. 



Pythia simulations, were approximated by the leading pion spectrum obtained in pion-proton 
collisions. Regeneration of A c -baryons was mimicked using an ordinary A as a beam particle in 
Pythia and extracting the spectrum of leading A's. The resulting muon and neutrino spectra 
are rather insensitive to these approximations, since they are slowly varying functions of kaon 
and heavier hadron regeneration Z-moments (they enter only through the combination Am in 
Eq. @) below). 

For the charm Z-moments one has to consider that the charm cross section need not scale 
with the target atomic number A in the same way as the total inelastic cross section. The 
ratio of the ^-dependence in our charm cross section (see section 5.2) and the inclusive A 2 / 3 - 
dependence (see section 2.3) gives a ^4 1//3 -dependence in the Z-moments for charm, which are 
included in our results (e.g. Fig. |j). 

The energy dependence of the hadron generation and regeneration Z-moments is shown in 
Fig.|3] for a constant spectral index 7 = 1.7 (dotted lines) and for a primary spectrum with 
a knee as in Eq. (|l|) (solid lines). For comparison, we also show the constant values given by 
Lipari Q under the assumptions of energy-independent inelastic cross sections and Feynman 
scaling of the meson spectra |34|, |35]| , i.e. with a Z- moment defined by Eq. (]24|). It is clear 
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Particle 


7T± K± 


K L 




D° 


A c D s 


J/1> 


e [GeV] 


115 850 


205 


3.74 • 10 Y 


9.47 • 10 Y 


2.57 • 10 8 9.33 • 10 Y 


8.8 • 10 15 



Table 2: Critical energies for hadrons traversing the atmosphere in the vertical direction. 



from the figures that variations of the Z-moments with energy are non-negligible, in particular 
when the changing slope (7) of the primary spectrum is included. On the other hand, this 
energy dependence of the Z-moments, and of the interaction lengths A, is mild with respect to 
the rapid decrease of the primary flux with increasing energy. For example, differentiation of 



Eq. (25) gives 

E d(f> N X E dAjy 



(7 + 1), (30) 



4>n dE An Aat dE 

where the first term comes from the energy dependence of Z and A and the second term from the 
primary spectrum. Numerically, the former is only 0.1, much smaller than the latter which is 
2.7 or 3.0. This demonstrates that the treatment of non-scaling effects as corrections embodied 
through energy-dependent Z-moments is justified irrespectively of the derivation leading to 
Eq.dl). 

We can now find approximate asymptotic solutions for the muon and neutrino fluxes. In 
the asymptotic regimes meson fluxes are well approximated by power laws, 4>m{E) oc E~&, with 
(3 = 7 in the low energy case and (5 = 7 + 1 in the high energy case. Hence the source terms in 
the lepton cascade equation (Eq.[l^)) can be rewritten as 

S(M - tY) = Z^+j-^-, (31) 

with decay spectrum- weighted moments defined byQ 

v (EmY P d M {E) dn(M^£Y;E M ,E) 

Zm ^ p+1 = J E dEM VW) -omem) d^ • (32) 

Integrating the lepton cascade equations over the line of sight, one then obtains the following 
expressions for the lepton fluxes deep in the atmosphere (X ^00): 

4™ = Z M -+i,t+i 1 Z ™ <Pn{E) (33) 

for leptons coming from low energy mesons and 

.high nr Z NM ln(A M /AAr) e M , 

h = Z M ^ a+2 T -^— i _ An/ ^ — ME) (34) 



for leptons coming from high energy mesons (with Em S> ttim)- I n Eq. (|34D, £m is a critical 
meson energy separating the low energy and the high energy regimes, i.e. where the meson 
dominantly decays or interacts, respectively. It depends on the atmospheric profile and in 



general on zenith angle. For the exponential atmospheric profile in sect. 2.2 and in the vertical 
direction, one has J|] 

m M c 2 h 

e M = • (35) 

CTM 



4 To keep the analogy with the (re)generation ^-moments, which include the multiplicity of the final state, we 
include the branching ratio Br(M — > IX) into the definition of the decay Z- moments. This differs from 
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0.00566 


D° — 




0.00839 


0.00636 


0.00354 


0.00283 


D s -» 




0.00744 


0.00550 


0.00292 


0.00228 


A c - 




0.00395 
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0.00298 
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0.0191 






0.0187 
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Table 3: Decay Z -moments Zm^£,/3+i for various decay channels. (For the decay J / 'ip — » \x 
a factor 2 is included to account for both [i~ and [i + .) 

In Table || we have collected the critical energies em used in this study. For not too large zenith 
angles 6^60°, Eq. @ leads to em oc 1/ cos# and (j)^ gh depends on the zenith angle. 

To resume, the lepton fluxes from meson M have the same spectral index of the primary 
flux, $ w oc ET*- 1 , and are independent of zenith angle at energies smaller than the meson 
critical energy em, while they are steeper by one power of energy and depend on zenith angle 
at energies above em, <j>i 9 oc I?~ 7 ~ 2 / cos#. We see that at the energies of interest to us, pions 
and kaons are above their critical energy and so generate 'conventional' muons and neutrinos, 
while charmed mesons are below their critical energy (~ 10 7 GeV) and give 'prompt' muons and 
neutrinos. 

The energy spectra of muons and neutrinos from decays of ultra-relativistic mesons take a 
simple scaling form || 

f E \ dE 

dn(M £Y;E M ,E) = F M J — — , (36) 

\EM J 

and the decay Z-moments are independent of energy 

Zm^I,P+i= [ dxx^F M _ e (x) , (37) 
Jo 

with x = E I 'Em- Approximate expressions for the functions Fm^i have been obtained for two 
and three body decay channels in ref. ||.|] Since there are many semi-leptonic decay channels 

5 Because of our convention, the functions Fm->i in ref. j^| should be multiplied by the branching ratio 
Br (M — ► £ X). 
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for charmed mesons and most of them have more than three particles in the final state, we prefer 
to generate all decay spectra within the Lund Monte Carlo. In Table || we list the values of the 
decay Z-moments for the spectral indices of interest in this study. 

Finally we join the low and high energy solutions with the interpolation 

\ - <p£° W 4>} lgh 4>n(E) Z NM Z M ^ t ^ + i /„o\ 

<pi™ + ^igh ~ i _ Znn 2_, ! + A M E/s M ' 



with 



Z M ^£ i7+ i 1 - A N /A M 

The fluxes of muons and neutrinos calculated according to Eq. ( |38[ ) using the previously- 
obtained energy-dependent Z-moments are plotted in Fig. || as dashed lines. It is satisfying to 
see that the cascade simulations and the approximate analytic solutions, which are conceptually 
rather different, give results that are quite close. Detailed comparison of corresponding fluxes 
shows good agreement both for conventional and prompt muons and neutrinos. The differences 
are typically less than 20% which is quite sufficient in this context. For the prompt leptons, 
this is below the uncertainty in our charm calculation (see section 5) and far smaller than the 
differences between the different models discussed in section [7|. 



5 The model for particle production 

A model for particle production is needed to specify the energy spectra of secondaries in cosmic 
ray collisions with atmospheric nuclei. As discussed in section |2.3| , collisions involving nuclei 
can be reduced to the simpler proton-nucleon collision. This applies in particular when only 
energetic particles are of interest, as in our case. 

The flux of conventional muons and neutrinos results from the decay of relatively long-lived 
particles, such as tt and K mesons. The production of such hadrons, containing only light quarks 
(u,d,s), is dominated by minimum bias proton-nucleon interactions (without large momentum 
transfers) and receives a small contribution from diffractive interactions. On the other hand, 
the prompt muons and neutrinos arise through decays of short-lived particles, i.e. dominantly 
charmed particles. Charm quarks are, due to their relatively large mass, usually considered 
to be produced in hard processes which can be described by perturbative QCD (pQCD). In 
the following, some relevant details of the models implemented in the Pythia and Jetset 
Monte Carlo programs will be discussed. The hypothetical non-perturbative intrinsic charm 
mechanism is discussed separately in section ||. 



5.1 Light particle production 

The production of light hadrons is dominantly through minimum bias hadron-hadron collisions. 
The strong interaction mechanism is here of a soft non-perturbative nature that cannot be cal- 



culated based on proper theory, but must be modelled. In the successful Lund model |3^] hadron 



production arise through the fragmentation of colour string fields between partons scattered in 



semi-soft QCD interactions [31]. The essentially one-dimensional colour field arising between 
separated colour charges is described by a one-dimensional flux tube whose dynamics is taken 
as that of a massless relativistic string. Quark-antiquark pairs are produced from the energy 
in the field through a quantum mechanical tunneling process. The string is thereby broken 
into smaller pieces with these new colour charges as endpoints and, as the process is iterated, 
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primary hadrons are formed. These obtain limited momenta transverse to the string (given by a 
Gaussian of a few hundred MeV width) but their longitudinal momentum may be large as it is 
given by a probability function in the fraction of the available energy-momentum in the string 
system taken by the hadron. The iterative and stochastic nature of the process is the basis for 
the implementation of the model in the Jetset program 31]. 



A non- negligible contribution to the inclusive cross section is given by diffractive interactions. 



These are also modeled in Pythia [31| using cross sections from a well functioning Regge-based 
approach and simulating the diffractively produced final state using an adaptation of the Lund 
string model. These diffractive events are included in our simulations and contribute rather less 
than 10% to the final results. 

5.2 Perturbative production of charm 

Charm production is strongly suppressed in hadronisation models that are commonly used in 
detailed comparisons with experimental data on multihadron production in various high energy 
collisions. The tunneling mechanism in the Lund model |^] gives a production probability of 
different quark flavours as uu : dd : ss : cc ~ 1 : 1 : 0.3 : 10 -11 , i.e. charm production in the 
hadronization phase can here be safely neglected. 

Charm quarks are instead considered to be produced in perturbative QCD processes in 
accordance with the relatively large charm quark mass. To leading order (LO) in the coupling 
constant, i.e. 0(a 2 ), these are the gluon-gluon fusion process gg — > cc and the quark-antiquark 
annihilation process qq — > cc as shown in Fig.|5]abc. The charm production cross section is 
calculated using the usual convolution of parton densities in the colliding hadrons and the 
hard parton level cross section a from pQCD, i.e. 

a = II J dx ^ di f^Q 2 ) h( x 2,Q 2 ) ^ (40) 

Here, Xi are the parton longitudinal momentum fractions in the hadrons and t is the Mandel- 
stam momentum transfer at the parton level. Q 2 is the factorization scale defining at what 
momentum transfer the parton densities are probed and also regulating the amount of pQCD 
scaling violations; we have used Q 2 = (?n^_ c + m^_-)/2, where = m 2 + p\. The charm 
quark mass introduces a threshold in the invariant mass of the parton level subsystem, i.e. 
s = X1X2S > 4m^c 4 . The dominating contribution to the cross section comes from the region 
close to this threshold, since dajds is a steeply falling distribution. It is therefore important 
to use QCD matrix elements with the charm quark mass explicitly included. The numerical 
value used is m c = 1.35 GeV/c 2 together with Aqcd = 0.25 GeV (in accordance with using the 
MRS G parton density parametrisation |3(| ) . 

Next-to-leading order (NLO), i.e. O(o%), cross sections for heavy flavour production in 
hadron collisions have been calculated in pQCD [37, Compared to the leading order re- 



sults there is an overall increase of the cross section of about a factor of two. This does not 
demonstrate a bad convergence of the perturbative series, since the main NLO contribution is 
associated with a process that does not appear in leading order charm production. This is the 
gluon scattering process gg — > gg, which has a much larger cross section than the leading order 
charm processes and is of a comparable magnitude when including the NLO correction g — > cc 
shown in Fig.[5]d. Since the NLO distributions of the charm quark transverse momentum and 
rapidity to a reasonable approximation have the same shape as the LO ones, we take the NLO 
results into account by rescaling the cross section with an overall factor K = 2. Still higher 
order corrections have not been calculated, but their effect should be significantly smaller since 



15 




a) 



b) 



c) 



d) 



Figure 5: Illustration of charm production in pQCD. The leading order processes (a,b,c) and an 
important next-to-leading order process (d). 



there should be no such additional process entering like for NLO. The factor K = 2 is indeed 
consistent when comparing the leading order cross sections with experimental results p9[ . 

Another estimate of higher order corrections can be obtained from charm production in 
the simulated parton cascades implemented in Pythia. These represent a leading logarithm 
approximation of multiple parton emission from the incoming and scattered partons in basic 
QCD 2^2 processes. Charm quark production arises here through the perturbative QCD 
gluon splitting process g — > cc, i.e. basically the same as in the NLO matrix element calculation. 
This approximate charm production to arbitrary order in pQCD gives a contribution of the same 
magnitude as the LO matrix elements, thus confirming the use of a renormalisation factor of 
K = 2. Furthermore, the energy spectra of the charm quarks from this higher order treatment 
are very similar to those from the LO calculation. If anything, they rather tend to be slightly 
softer than the LO distribution Q. Since it is the hardest part of the spectrum that gives 
the largest contribution to the high energy neutrino and muon spectra, one may conclude that 
taking the higher order corrections into account through a global i^T-factor renormalisation and 
keeping the leading-order shape of the charm quark energy spectra, as we have done, is sufficient 
for the precision needed in this study. 

The charm production cross section depends on the input for the parton density functions 
fi(%i, Q 2 ) in Eq. (|40|) . With charm production being dominantly close to threshold s = X1X2S > 
4m;? c 4 , the typical initial momentum fractions xi will decrease with increasing collision energy s. 
This is demonstrated in Fig.||, which shows the distribution of initial parton momentum fractions 
in charm production at different energies. At the highest energies, the parton densities are probed 
down to x ~ 10 -5 or even below. The recent data from the ep collider HERA [40, 41 1 show a 
significant increase at small x, xf(x) ~ x~ a and constrain the parton densities down to x ~ 10 . 
These data, together with other data from previous deep inelastic scattering experiment as well 
as other processes, have been used in the parametrization MRS G [^] of parton densities. The 
resulting small- x behaviour is given by the power o = 0.07 for sea quarks and a = 0.30 for gluons. 
Since MRSG is the most recent parametrisation, using essentially all relevant experimental 
data, we use this as our standard choice. To investigate |l the dependence on the choice of 



parton density parametrisations, we also applied the MRS Dq [42| with the small- x behavior 
x f(x) ~ const, which before the HERA data was an acceptable parametrisation. The effect 
on the total charm production cross section from the choice of parton density parametrisation 
is illustrated in Fig.[7| with curves resulting from these parametrisations. At high energy there 
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Figure 6: Distribution of momentum fraction x for the initial partons entering the QCD charm 
production subprocesses. The curves represent three different beam particle energies: 10 3 GeV 
(solid), 10 6 GeV (dashed) and 10 9 GeV (dotted). 



is a large dependence on the choice of parton density functions. The difference between the G 
and the Dq parametrisations should however not be taken as a theoretical uncertainty. First of 
all the Dq parametrisation is known to be significantly below the small- x HERA data and gives 
therefore a significant underestimate at large energies. Secondly, the naive extrapolation of the 
G parametrisation below the measured region x^lO -4 at rather small Q 2 (~ m 2 ) leads to an 
overestimate. A flatter dependence like x~ e with e ~ 0.08 as x — > can be motivated ( fS^ and 
references therein) based on a connection to the high energy behaviour of cross sections in the 
Regge framework. The implementation of this approach in Pythia makes a smooth transition 
to this dependence such that the parton densities are substantially lowered for x^lO -4 leading 
to a substantial reduction of the charm cross section at large energies, as given by the solid curve 
in Fig.0. 

Using this procedure we get a quite decent agreement between experimental charm produc- 
tion data and the Pythia simulation results over a wide range of energies (Fig. [7]). A few 
comments on the data in Fig. ^ are here in order. A given experiment is only sensitive to some 
channels and a limited kinematical region. The total charm cross section is therefore obtained 
by a rescaling with charm decay branching ratios and by using assumed shapes of the xf dis- 
tributions to extrapolate to unmeasured regions. In particular, corrections to points 1,2,6 and 
7 are small while they are large for point 9 and 13. The bands 8,10,11 and 12 illustrates the 
uncertainty in these experiments due to this extrapolation. In band 8 the uncertainty includes 
a scaling for including D^-mesons (taken from ^9|, 51]). Data-band 11 is based on D + D iden- 
tification, 12 on AjD and 10 on A+ identification. Furthermore, points 3, 4 and 5 are from 
beam dump experiments on heavy nuclear targets without direct charm identification and have 
an additional uncertainty from the scaling with nuclear number. In point 4 the scaling A 
has been assumed, which we have rescaled in 4' to a ^-dependence in order to be consistent 
with the other beam dump experiments and with our model. Data points 2 and 6 come from 
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Figure 7: Energy dependence of charm production cross section in pp(pp) . The experimental data 
points are 2:jm, 3:g§j, 4,4'^, 5:^, 6:g$, 7:j$J, 8:^J 9:^, 10:^, \5$, 

12:^5b\, |5^/ and 13:^5%j. The solid line is the result of our model. The dotted lines result from 
a naive application of the MRS parametrisations G and Dq of parton densities. The dashed 

which is 10% of the 



lines represent earlier models, C from \ldjl , V from t V7j and Zl from 
total pp cross section a tot . The curves IC1 and IC2 are based on the intrinsic charm hypothesis 
(section^. Detailed discussions are in sections 5.i and^. 



pp interactions with explicit charm particle identification. Although these issues leave some 
uncertainty for each individual result, the combination of all data should give a trustworthy 
knowledge on the charm cross section and its energy dependence. 

As indicated, there is an uncertainty related to the dependence of the charm production 
cross section on the nuclear mass A. Since the discussed pQCD charm production involves 
hard scattering processes with a small cross section, one may argue that it does not only take 
place with nucleons at the nuclear surface but also with nucleons in the interior. Consequently, 
the power a c in a(pA — > cc) = A ac o~(pN — > cc) may be larger than 2/3 and rather be 1, 
corresponding to interactions in the whole nucleus. In fact, experiments give values of a c that 



are mainly close to 1 (see e.g. [63] and references therein), which is the value we have adopted 
in our calculations. This is also in agreement with a number of recent experiments with both 
pion and proton beams [^l], 58, 59 1. Using a c = 2/3 instead would lead to a reduction of the 
normalisation for the prompt muon and neutrino fluxes by the factor A 2 ^ 3 /A = 14. 5 -1 / 3 ~ 0.4. 
This is at most an upper bound for the uncertainty of the ^-dependence. 



5.3 Semileptonic decays of charmed hadrons 

The branching ratios for semileptonic charm decays used are BR{D^ —* e/fj,) = 17%, BR(D° 

e/fi) = 8%, BR(Df -» e/n) = 8%, BR(A+ -» e/y) = 4.5%. These values agree with experi- 
mental measurements as compiled by the 'Particle Data Group' [|33f] . To distribute momenta in 
the charm hadron decay H — > iv^h, the Monte Carlo [31] uses the weak V — A matrix element 



\M\ 2 = (p H Pi)(PuPh) 



(41) 
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in terms of the involved four-vectors and neglecting decay product masses. H is the charm 
hadron and h is the final hadron in a three-body decay and generalized to represent the final 
state hadron system in case of more final hadrons. The multiplicity and momenta within the 
/Vsystem is phenomenologically modelled |H| and tested against data. These details, however, 
concern the final hadrons and are not important for our purposes, because we only use the final 
leptons, which should be well accounted for through the weak decay matrix element (O) and 
the known branching ratios. 



6 Possible non-perturbative origin of charm 

Although most of the charm production data from accelerator experiments can be reasonably 
well understood from pQCD calculations, the uncertainty in the data and the calculations cannot 
exclude some smaller non-perturbative contribution. Charm production in pQCD is theoreti- 
cally well defined and only has some limited numerical uncertainty due to parameter values 
and NLO corrections which, however, can be examined and controlled as discussed above. In 
contrast, non-perturbative charm production is not theoretically well defined due to the general 
problems of non-perturbative QCD. In particular, the absolute normalisation has to be taken 
from comparison with data. Models for non-perturbative charm production exist and some 
have been used in the context of atmospheric muon and neutrino fluxes, as discussed in sec- 
tion 0. Here, we will investigate the consequences of the hypothesis of intrinsic charm quarks in 
the nucleon wave function poL Although being far from established, this idea has theoretical 



motivation fl60| , 61] and some experimental data can be interpreted as giving evidence [62-67]. 
It is therefore a serious model for a possible additional mechanism for charm production with 
a non-perturbative origin. The results presented in this section are based on a more general 
study HI of intrinsic charm in high energy collisions, to which we refer for more details of our 
implementation in an explicit model. 

The hypothesis of intrinsic charm (IC) amounts to assuming the existence of a cc-pair as a 
non-perturbative component in the bound state nucleon |6(|. This means that the Fock-state 
decomposition of, e.g., the proton wave function, \p) = A\uud) + B\uudcc) + contains a small, 
but finite, probability B 2 for such an intrinsic quark-antiquark pair. This should be viewed as a 
quantum fluctuation of the proton state. The normalization of the heavy quark Fock component 
is the key unknown, although it should decrease as I/tuq. Originally, a 1% probability for charm 



was assumed, but later investigations, e.g. [64, ^5[, indicate a smaller but non- vanishing level. 

Viewed in an infinite momentum frame, all non-perturbative and thereby long-lived compo- 
nents must move with essentially the same velocity in order that the proton can 'stay together' 
for an appreciable time. The larger mass of the charmed quarks then implies that they take 
a larger fraction of the proton momentum. This can be quantified by applying old-fashioned 
perturbation theory to obtain the momentum distribution |6 



P(p — > uudcc) cx 



5 ™2 



771. i 



1=1 Xi 



(42) 



in terms of the fractional momenta X{ of the five partons i in the uudcc state. Neglecting the 
transverse masses of the light quarks in comparison to the charm quark mass results in the 
momentum distribution 



r2 ™2 



P(xi,x 2 ,x 3 ,x c ,Xc) oc - — c c _x 2 £(l -xi-x 2 -x 3 -x c - xi) (43) 
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which favour large charm quark momenta as anticipated. In fact, one obtains (x c ) = 2/7 by 
integrating out the light quark degrees of freedom x%. 

A proton with such an intrinsic cc quantum fluctuation can then interact with another 
hadron such that charmed particles are realised. A hard interaction with such a charm quark 
is certainly possible, but the cross section is then suppressed both by the small probability of 
the fluctuation itself and by the smallness of the perturbative QCD interaction. The charm 
quarks may, however, also be put on shell through non-perturbative interactions that are not 
strongly suppressed fl6l| . This may lead to a rate that is large enough to be of potential interest. 
To estimate these non-perturbative interactions we have constructed a model p8| based on 



refs. J60|, |62 



The formation of charm hadrons can occur through the following mechanisms. The charm 
(anti)quark can hadronise into a Z}-meson as described by a normal hadronisation function, 
similar to those successfully used in e + e~ — > cc. Alternatively, the (anti)charm quark can 
coalesce with another quark or diquark from the \uudcc) state to form a hadron. Following 
fl62"| , l&f l we use the recombination probabilities 50% to form a L>-meson and 30% for a A c , 
in which cases the remaining core quark is assumed to hadronise separately from the proton 
remnant. The probability to directly form a J/ip {i.e. the cc pair is combined) is taken to 
be 1%. The momentum of the hadron formed through coalescence is taken as the sum of the 
corresponding Xj's, e.g. x\ c = x c + x u + xj. The momentum distribution is then obtained by 
folding Eq. (|^) with the proper 5 function, e.g. S(xa c — x c — x u — Xd), and integrating out all 
extra degrees of freedom. The core quarks that does not coalesce with spectator partons are 
hadronised to L>-mesons with a normal hadronisation function. Since such a function is quite 
hard, we here approximate it with a ^-function to let the charm hadron take the whole charm 
quark momentum given by Eq. (^), which is consistent with low-pt charmed hadroproduction 
data [ ]67| . 

The shapes of the ^-distributions for the charmed particles are thereby given (see |6q| ). 
They are quite hard, in fact harder than those for charm from pQCD, and therefore have 
the potential to contribute effectively at high energies. As mentioned, the main uncertainty 
in the intrinsic charm model is the absolute normalization of the cross section and its energy 



dependence. The magnitude of the cross section has been estimated [63] from data at relatively 
low energies (E p = 200 — 400 GeV). Since the process is basically a soft non-perturbative process 
it may be reasonable to assume that its energy dependence is the same as that for normal inelastic 



scattering [61, p8fl . We therefore take as our first case 

ICl: a IC {s) = 3 • W- 5 a mel (s) (44) 

with normalisation from |33| and shown as curve ICl in Fig.|7[ Alternatively, one might argue 
that there is a stronger energy dependence related to some threshold behavior for putting the 
charm quarks on their mass shell. We make a very crude model for this by taking the intrinsic 
charm cross section to be a constant fraction of the pQCD charm cross section 

IC2: <ric(s) = 0.1 a pQCD (s) (45) 

as shown by curve IC2 in Fig. 0. This is similar to the low energy (200 — 800 GeV) treate- 
ment in JS^ [. The normalisation is here fixed to be the same as ICl at the low energy where 
evidence is claimed for intrinsic charm [^]. There is, however, some indication against such an 
increased cross section, as in IC2, since no evidence for J ftp from intrinsic charm was found in 
an experiment |39| at a somewhat higher energy (800 GeV proton beam). 

The dependence on nuclear mass number should (in both cases) essentially be A 2 / 3 reflecting 
the soft nature of the hadron-hadron interaction. Note, that the intrinsic charm quarks can be 
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Figure 8: Energy- dependence of production Z^-moments, Eq. $2!\), for incoming particle 
k = proton producing charmed hadron h = Z)°,£) ± ,A C and J/i/j through the intrinsic charm 
mechanism with the two assumed energy dependences IC1 and IC2. 



released through interactions with the other parts of the \uudcc) state, as demonstrated clearly in 
case of J/ip production after the remaining proton state has been 'stripped off' in the interaction 

The intrinsic charm model provides simple and scaling x^-spectra for the charmed particles 
which makes the analytic method suitable for calculating the fluxes of muons and neutrinos 
from their decays. The charm production Z-moments are calculated according to Eq. (p9|) by 
numerical integration and shown in Fig. ||. The mild energy dependence of the IC1 case gives 
essentially constant Z-moments, except for a marked step due to the change of the slope 7 in 
the primary spectrum. The step is more pronounced here compared to charm from pQCD (see 
Fig.Q) where the ^-distribution is both energy-dependent and extending to smaller values. 
The stronger energy dependence in IC2 is reflected in the corresponding Z-moments, where the 
strong variation at energies below ~ 10 5 GeV may cast some doubt on the reliability of the 
analytic method in this region. 

The lepton fluxes are then obtained by using Eq. ( |38[ ) and the regeneration and decay Z- 
moments calculated in section 4. The results are displayed in Fig.^ and compared to those 
from our pQCD calculation. The milder, and more conservative, energy dependence (IC1) of 
the intrinsic charm cross section gives a result which is only a small (~ 10%) correction to the 
pQCD result, except at super-high energies where the rate is extremely small and not measurable 
in a foreseeable future. Note that the J/ift contribution is here becoming important, since the 
J/ip flux is not attenuated through interactions due to the high critical energy (Table Q). This 
raises the question of how well the normalisation of the otherwise small J/ip contribution is 
known. With the strong energy dependence assumed in IC2, the intrinsic charm result exceeds 
the pQCD one already at lepton energies around 10 4 GeV. Although the energy dependence of 
this IC2-model is, as mentioned, rather ad hoc and may be disfavoured by data, it illustrates 
the large theoretical uncertainty associated with the intrinsic charm model when extrapolated 
to the high energies of cosmic ray interactions. 



21 




IC2 



, prom pt 




8 10 2 4 

log 10 (E/GeV) 



' 



Figure 9: Fluxes (cf. fig. 3) of muons and neutrinos from the intrinsic charm model with the 
two assumed energy dependences IC1 and IC2; for IC1 the contributions from different charm 
particles are given. Shown for comparison are our standard results (full curves) for conventional 



and prompt fluxes as given by the parametrisation Eq. and extrapolated to higher energies 
(dashed). 



7 Comparison with previous model calculations 

In sections 3 and 4 we have obtained atmospheric muon and neutrino fluxes with two different 
methods: via a Monte Carlo simulation of the hadronic cascade and via approximate analytical 
expressions with energy-dependent Z-moments. We were satisfied that the two methods gave 
consistent results. Here we want to compare our results with those obtained with different 
models for particle interactions in the atmosphere. In particular, we focus on the prompt muon 
and neutrino fluxes arising from different charm production mechanisms. 

Earlier calculations of the conventional muon and neutrino fluxes jj], [?], ||, ||] agree well with 



our results as shown in Fig. 1C . The conventional muon flux from Gaisser O is shown as a dashed 



line in Fig.QOk as far in energy as it is applicable. Also shown as dashed lines in Fig. 10 be are the 



conventional neutrino fluxes from Volkova M. In contrast to these models, ours does not obey 



Feynman scaling. The scaling violations are apparent in Fig.|ll|a and from the deviation from 
a constant value of the production Z-moments (Fig. Q). Nevertheless, they are small enough 
that when folding everything together (initial spectrum, cascade interactions and decays) the 
resulting conventional muon and neutrino fluxes agree well with those models. We have thus 
confirmed previous results by an independent calculation based on a new approach using Monte 
Carlo simulations to more fully take into account the atmospheric cascade interactions producing 
secondary particles decaying into muons and neutrinos. 

Concerning previous estimates of the flux of prompt muons and neutrinos, there are variations 
between different model calculations of up to a few orders of magnitude, as illustrated in Fig. |n]. 
One should note that prompt fluxes are direction independent up to the charm particle critical 
energy ~ 10 7 GeV (Table ^) and therefore directly comparable independently of whether the 
horizontal or vertical direction has been considered in these estimates. Furthermore, due to the 
charmed particle decay kinematics and the same branching ratios for the semi-leptonic decays 
into electrons and muons, the prompt muon and neutrino fluxes are essentially the same (cf. the 
decay Z-factors in Table 0). Therefore, the curves for prompt muons in Fig.[To|a are also taken 
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Figure 10: Prompt and conventional muon and neutrino fluxes from our cascade simulation 
(solid curves) compared with earlier model calculations as discussed in the text. 



to represent the prompt neutrinos in Fig. [To|bc (except for our own curves, which are calculated 
separately) . The comparison in Fig. |l0| shows that previous results are in general substantially 
larger than our results based on pQCD. Even the highest flux from our extreme version IC2 
of the intrinsic charm hypothesis (Fig. ^) is lower than most previous calculations. These large 
differences are due to different models for charm production, both regarding the magnitude 
and the energy dependence of the cross section and the distribution in longitudinal momentum 
fraction xf of the charmed particles. 



The curves labeled V in Fig. 10 are from the calculation by Volkova et al. 17], applying 
the so-called 'quark-gluon string model' (QGSM) |Fo| (not to be confused with the Lund string 
model f32|]). It uses a parametrised energy dependence of the charm cross section, curve labeled 
V in Fig.[7|, normalized to early experimental data which are substantially above more recent 
measurements. The charm particle energy spectrum is assumed to obey Feynman scaling and 
has the form 

dN/dx F ~ (1 -x F ) a (46) 
with old = 5 and a\ c = 0.4 for D-mesons and A c -baryons, respectively. 



The curves marked Zl in Fig. 10 are from Zas et al. |2(| and illustrates an extreme model 
where the charm cross section is simply taken as 10% of the total inelastic cross section 
(cf. Fig.[7|). This is substantially higher than all charm data as shown in Fig. ^. This model 
uses the scaling x^-distribution of Eq. (H) with = 3 and a^ c = 1- Castagnoli et al. [jD| 
obtained the result marked C in Fig. |l^ using a parametrised energy dependent charm cross sec- 
tion shown by curve C in Fig.f?] based on some early data (band marked 11) that are higher than 
later measurements. Again, the differential spectra are of the form Eq. ( |46| ) using an = 5 and 
a\ c = 0.4. The curves marked Z2, from Zas et al. |po"| , correspond to charm quark production 
calculated with leading order pQCD matrix elements using relatively hard parton distributions. 
This spectrum would be softened by taking hadronisation into account and thereby become even 
closer to our result, as expected since they are based on the same pQCD processes. 

The first important difference between our model and previous ones lies in the magnitude 
and energy dependence of the charm production cross section. As demonstrated in Fig.|?] our 
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Figure 11: Distribution in fractional energy xe of produced mesons: (a) tt and mesons, (b) 
D mesons. Dashed, dotted and dashed-dotted lines are obtained from simulations with Pythia at 
proton beam energies of 10 3 GeV, 10 6 GeV and 10 9 GeV, respectively. The solid line represents 



Feynman scaling using Eq. (46) with old = 5. (D-meson curves are normalized to unit area.) 



model reproduces available data on charm production cross sections, but the other models do 
not. In some cases one may have been mislead in the construction of the models by the early 
charm measurements that turned out to be substantially higher than the measurements done 
later. 

Another important reason for our lower flux is the strong breaking of Feynman scaling 



as demonstrated in Fig. lib. The xf distributions for D-mesons produced at three different 
energies in our model are here compared with the energy-independent distribution in Eq. ( |4"6| ) 
with old = 5. The Feynman scale breaking in our model arises in the perturbative charm 
quark production, but is also influenced by the hadronization model. As discussed in section 
5.2, charm quark production is dominant close to threshold, s = x±X2S > Am^c 4 . This effect 
does not disappear with increasing energy, but is rather enhanced with parton densities that 
increase at small x. This leads to a scale breaking with the charm quark xf distribution in the 
symmetric nucleon-nucleon cms becoming softer around xf = with increasing cms energy. In 
the Lund hadronization model, the charm quark is connected by a colour string to a spectator 
parton. In the hadronization of this string, the produced charm hadron may obtain a larger 
longitudinal momentum than the charm quark, due to the momentum contribution from the 
parton with which it is joined. The string may even have so small invariant mass that it directly 
produces a charmed hadron, i.e. the charm quark effectively coalesces with a spectator parton 
into a charmed hadron. This latter process naturally happens particularly at small overall 
cms energies. Since these forward- 'pulling' hadronization effects become less important with 
increasing collision energy, they also contribute to the Feynman scale breaking. 

The effect on the prompt muon flux from this Feynman scale breaking is shown in Fig. |l2| 
Our normal result is here compared with the results from a modification of our model, where in 
each event the charmed particles D and A c are redistributed according to the scaling distribution 
in Eq. ( |46| ) with old = 5 and ot\ c = 0.4. Clearly, the Feynman scale breaking softens the spectrum 
considerably. (In comparing these curves one should note that there is no conserved integral of 
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the E 3 -weighted flux.) To examine the effect of the energy dependence of the overall charm cross 
section, we have, in addition to using this scaling a^-distribution, renormalized our simulated 
charm events to mimic the cross section in the model by Volkova et al. [17] mentioned above 
and shown by curve V in Fig. ^. Since this cross section is larger than in our model at energies 
below ~ 10 6 GeV and smaller above, this change flattens the muon spectrum in Fig. 12. With 



these two changes in the spirit of ref. [17|, we obtain the same shape of the prompt muon flux as 
in 1 17] but with a lower overall normalization. The resulting spectra are, however, in reasonable 
agreement with the calculation of Castagnoli et al. [15] using a similar approach as ref. [17]. A 
more recent calculation based on the QGSM by Gonzalez-Garcia et al. [21] gives fluxes that are 
comparable to the fluxes predicted by Castagnoli et al. [15|. 

The calculation by Bugaev et al. |lq 1 resulted in overall prompt fluxes slightly lower than 
in ref. [17|. They considered Feynman scaling violations in charm production through a phe- 
nomenological equation and obtained higher fluxes in the non-scaling case than in the scaling 
case, i.e. opposite to the effect we find and have just described. However, their way of introduc- 
ing the energy dependence in the xp distribution does not preserve the overall normalization, 
i.e. the integral of the ^-distribution. This means that, in comparison with our model, there 
is not the same clear separation between the overall charm cross section normalization and the 
charm particle xp distribution. 

Since most of the earlier calculations are based on non-perturbative charm production mech- 
anisms a comparison with our intrinsic charm model in section || is of interest, i.e. comparing 
Figs.p and 1C. Intrinsic charm should be considered as a process in addition to the standard 
pQCD one, but the sum of their resulting fluxes is still lower than, e.g., Volkova et al. [17]. To 
get a similarly high flux the rate of intrinsic charm must be increased substantially, about a 
factor 1000 for IC1 and 10 for IC2. Such rates are incompatible with the experimental limits on 
intrinsic charm and with measured inclusive charm cross sections (Fig. [7]) and x^-distributions. 

This discussion has demonstrated significant effects on the high energy prompt muons and 
neutrinos depending on the assumptions made in the charm production model employed. The 




Figure 12: The prompt muon flux in our model (solid line) and after re- distributing the generated 
charm particles to obey Feynman scaling using Eq. (dashed line) and then also renormalizing 
the cross section to that of curve V in Fig.^ from [F\ j (dotted curve). 
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models in ref. [15, 17 ] give cross sections above more recent charm production data (Fig. 0) and 



apply simple Feynman scaling ^^-distributions. The model used by us, on the other hand, gives a 
fair description of measured charm production cross sections (Fig.[7|) and applies well-motivated 
charm particle momentum distributions with significant Feynman scaling violations. 

8 Conclusions and outlook 

We have studied the production of neutrinos and muons in the atmosphere by collisions of cos- 
mic rays with air nuclei, paying special attention to muons and neutrinos coming from decays of 
charmed particles (prompt fluxes). Two methods have been used to calculate the fluxes: a Monte 
Carlo simulation of the hadronic cascade in the atmosphere and an interpolation of asymptotic 
solutions to the transport equations. In both methods the Pythia Monte Carlo program has 
been used to simulate the primary collision and following cascade interactions. Results for the 
two methods are consistent. They agree with previous calculations of the conventional fluxes 
from decays of pions and kaons, but give substantially lower prompt components. This is due to 
different models for charm production, both regarding the energy dependence of the cross section 
and the longitudinal momentum distribution of the charmed particles. Whereas previous models 
give charm production cross sections above collected recent data and apply Feynman scaling for 
the longitudinal momentum distributions, our model gives a fair description of measured charm 
production cross sections and applies well-motivated charm particle momentum distributions 
with significant Feynman scaling violations. There is still some uncertainty (as discussed in 
section 5.2) when extrapolating this charm cross section calculation to the very high energies 
(10 9 GeV) needed for this study. Here, one cannot at present exclude a non-negligible contri- 
bution from some unconventional non-perturbative production mechanism. We investigate one 
such mechanism, namely that of intrinsic charm which has some theoretical motivation and 
some indications from data. This mechanism is likely to give a very small contribution to the 
total charm cross section, but the poorly known normalisation and energy dependence prevents 
a reliable prediction. Although disfavoured, a contribution ~ 10% of the pQCD charm cross 
section is presently not excluded, which through the harder momentum spectrum would lead to 
a dominant contribution of leptons at very high energies. 

We find that prompt muons and muon-neutrinos overcome the conventional fluxes at an 
energy of 10 6 GeV, which is substantially higher than in some earlier estimates. According 
to our results, it will therefore be harder to use measurements of the prompt atmospheric 
fluxes to estimate the total charm production cross section at high energy. The situation is 
slightly different in the case of the electron-neutrinos, for which prompt fluxes dominate above 
10 5 GeV. The electron-neutrino flux is, however, experimentally more difficult to measure and it 
is therefore a challenge to obtain the data needed to derive the charm production cross section. 

On the positive side, the lower atmospheric neutrino fluxes we predict are a less severe 
background to measurements of neutrinos from astrophysical sources (for a review on these 



see |71|). To illustrate this, we show in Fig. 13 the vertical fluxes of conventional and prompt 



atmospheric muon-neutrinos calculated by us (solid lines) together with expected neutrino fluxes 
from such sources. Cosmic ray interactions with the interstellar medium produce neutrino fluxes 
through processes similar to the atmospheric case and we show results derived from |72| in the 
direction of the galactic center (dashed upper curve) and orthogonal to the galactic plane (dashed 
lower curve). Two estimates of diffuse neutrino fluxes from active galactic nuclei are also shown 
(dotted line from ref. [fnj and dash-dotted line from ref. fT4|). At high energies (^10 5 GeV), all 
of these fluxes are in excess of our predicted atmospheric neutrino background. This provides 
interesting prospects for large scale neutrino telescopes to detect high energy neutrinos from 
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Figure 13: Vertical fluxes of conventional and prompt atmospheric muon-neutrinos from our 
simulation (solid lines) compared to some astrophysical sources: the flux from cosmic ray inter- 
actions with the interstellar medium as derived from ( dashed upper curve: in the direction 
of the galactic center; dashed lower curve: orthogonal to the galactic plane) and the estimated 
diffuse fluxes from active galactic nuclei (dotted line Jffij , dash-dotted line ftf$). 



cosmic sources. 
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